Correlation effects in bistability at the nanoscale: steady state and beyond 
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The possibility of finding multistability in the density and current of an interacting nanoscale 
junction coupled to semi-infinite leads is studied at various levels of approximation. The system is 
driven out of equilibrium by an external bias and the non-equilibrium properties are determined by 
real-time propagation using both time-dependent density functional theory (TDDFT) and many- 
body perturbation theory (MBPT). In TDDFT the exchange-correlation effects are described within 
a recently proposed adiabatic local density approximation (ALDA). In MBPT the electron-electron 
interaction is incorporated in a many-body self-energy which is then approximated at the Hartree- 
Fock (HF), second-Born (2B) and GW level. Assuming the existence of a steady-state and solving 
directly the steady-state equations we find multiple solutions in the HF approximation and within 
the ALDA. In these cases we investigate if and how these solutions can be reached through time 
evolution and how to reversibly switch between them. We further show that for the same cases the 
inclusion of dynamical correlation effects suppresses bistability. 

PACS numbers: 72. 10.Bg,71.10.-w,31.15.xm,31. 15.ee 



I. INTRODUCTION 

The phenomenon of bistability has been the subject 
of several studies in the field of quantum transport. In 
their seminal paper Goldman et al^ reported the obser- 
vation of bistability in the I-V curve of double-barrier 
resonant tunneling (DBRT) structures, thus stimulat ing 
many theoreticaP"^ and experimental investigational 
on the subject. The bistability is a non-linear effect in- 
duced by the electrostatic charge build-up in the quan- 
tum well and occurs in the bias window of negative dif- 
ferential resistance^ From the theoretical point of view, 
various techniques have been used to capture this phe- 
nomenon, ranging from a crude estimate of the charge 
build-up^ to self-consistent calculations at a mean field- 

leveiPfflmni 

With the increasing interest in transport through 
nanoscale devices, in particular using molecules as a pos- 
sible component of future electronic circuits, the study 
of intrinsic bistability in nanoelectronics has gained new 
attention. The possibility of finding molecular devices 
equivalent to conventional non-linear devices, such as 
diodes and transistors, is indeed an attractive perspec- 
tive. There have already been some successful attempts 
along these lines. Molecular devices with large on-off 
current ratio and large negative differential resistance, 
behaving simil arly to mesoscopic DBRT structures, have 
been reported! 11 * 12 ! So far, the great majority of bistabil- 



ity studies have been limited to the steady-state regime 
and performed within the framework of the Landauer 
formalism combined with static density functional the- 
ory (DFT) . At the Hartree level, bist abilit y was reported 
for a double quantum dot structure! 13114 1 In the context 
of time-dependent (TD) DFT, the inclusion of mem- 
ory effects beyond the adiabatic approximation is not 
straightforward and the development of accurate func- 
tionals to be used in numerical calculations is still un- 
der way. A promising and timely, even though compu- 
tationally demanding, alternative is the solution of the 
Kadanoff-Baym (KB) equations^" 2 ^ using self-energies 
from many-body perturbation theory (MBPT). The ad- 
vantage of MBPT over TDDFT is the inclusion of dy- 
namical exchange-correlation (XC) effects, i.e., effects 
arising from a frequency-dependent self-energy, in a more 
systematic way through the selection of suitable Feynman 
diagrams. Thus, MBPT provides an important tool to 
go beyond the commonly used adiabatic approximations 
and to quantify the importance of dynamical XC effects. 

The fundamental issue which we address in this pa- 
per is whether the bistability phenomenon found in static 
DFT, Hartree and Hartree-Fock (HF) approximation sur- 
vives when dynamical XC effects are taken into account. 
In contrast to DFT and mean-field approximations the 
steady-state equations of MBPT do not form a closed 
set of equations for the density only. This difficulty 
renders the search for the bistability regime in MBPT 
computationally very costly. To overcome the problem 
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we implement a time-dependent strategyf^HMSIHSll w e 
first solve the steady-state equations of DFT and mean- 
field theory to determine the parameter range for testa- 
bility. Then we go beyond the current state-of-the-art 
and provide a TD description of the bistability phe- 
nomenon in adiabatic TDDFT^SHU anc j TD mean field 
theory. We show how to switch between different sta- 
ble states by means of ultrafast gate voltages or exter- 
nal biases (driving fields). The possibility of reversibly 
switching between different stable steady-states is an as- 
pect that has remained largely unexplored Knowing 
how to steer the electron dynamics in real time we use the 
same driving fields in correlated MBPT simulations. The 
calculations are performed with the fully self-consistent 
second-Born (2B) and GW approximations which have 
recently been shown 3 -^ to agree with numerically exact 
TD-DMRG results. 35 In all cases studied here where adi- 
abatic DFT and HF theory predict bistability dynamical 
XC effects destroy the phenomenon. 

The paper is organized as follows: In Section [IT] we 
introduce the model used in our study. In Section III 
we introduce the TDDFT approach to transport, with 
a particular emphasis on the real-time propagation to 
study time-dependent transport phenomena. Assuming 
that the system reaches a steady state, the steady-state 
density can be calculated without explicitly propagating 
in time, as explained in Section |III B| In this Section 
a separate (fixed-point) analysis is given to determine 
whether or not a solution is stable. As an alternative 
real-time approach to transport, we introduce in Section 
|IV| the MBPT approach based on the solution of the KB 
equations. In Section [VJ we present the results of our 
numerical simulations for a certain set of parameters for 
which multiple solutions are observed in DFT. Finally, 
conclusions are drawn in Section IVIl 



II. THE MODEL 

In this work we study multistability in quantum trans- 
port, i.e., we study the question if, for a given set of pa- 
rameters, a biased system can have more than one steady 
state. Apart from proposing a way to classify the differ- 
ent steady states as stable or unstable, we also investigate 
by explicit time evolution the possibility to reversibly 
switch between different steady states by application of 
a suitably chosen, time-dependent perturbation. 

These questions are addressed in model systems. We 
consider a nanoscale device consisting of a few interact- 
ing sites contacted to two non-interacting tight-binding 
leads. Initially, the contacted system is in equilibrium 
at a given temperature and chemical potential. At time 
to = we switch on a bias in the leads and follow the 
time evolution of the perturbed system. 

The Hamiltonian of the system consists of three differ- 
ent parts 



where the Hamiltonian Hc{t) of the central device de- 
scribes a chain of Nq sites with a Hubbard-type on-site 
electron-electron interaction: 

N c No-l 

H c (t) = J2e?(t)dld ia - {VcQJi+io + H.c.) 
1=1 i=i 



(2) 
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Here, d\ a (d,v) denote creation (annihilation) opera- 
tors for electrons with spin a at site i. The ef(t) 
are on-site energies which may consist of an arbitrary 
time-dependent part, denoted as V g i(t), plus a time- 
independent part, ef . The nearest-neighbor hopping in 
the chain is Vc and U is the on-site Hubbard interaction. 

The non-interacting left (L) and right (R) leads, 
a = L, R, are described by one-dimensional semi-infinite 
chains with Hamiltonian 

oo 

7 
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E (vAacA+i™ + H.c) , (3) 



i=i 



H(t) = H c (t)+ E H a (t) + H 7 



(1) 



a=L,R 



with creation (annihilation) operators c\ aa (ci CTQ ) for elec- 
trons with spin a at site i in the lead a. The on-site 
energies e a and the hopping matrix elements V a are in- 
dependent of time and site index while W a (t) describes 
a time-dependent, site-independent bias applied to the 
lead a. 

Finally, the tunneling Hamiltonian which couples the 
leads to the device is given by 

H T = ~E \ Vlinkd[ a ClaL + Vimkdjv CC rCi (T ,R + H., 

(7 

(4) 

where we have adopted the convention that the site in the 
lead a connected to the device is labeled by the site index 
1 and Vii n k is the hopping between the central region and 
the leads. 



III. MULTISTABILITY IN TIME-DEPENDENT 
DENSITY FUNCTIONAL THEORY 

In this section we introduce the TDDFT approach®] 
which we will use to investigate multistability in our 
model system. In TDDFT, the complicated problem 
of interacting electrons is mapped, in principle exactly, 
on the much simpler problem of non-interacting Kohn- 
Sham (KS) electrons moving in an effective local poten- 
tial, thus providing a natural way to acco unt fo r correla- 
tion effects in both leads and device regionP3Hl Since the 
method only involves the propagation of single-particle 
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wavef unctions, it promises for a computationally efficient 
way to study time-dependent phenomena in quantum 
transport.^ The real-time TDDFT approach to trans- 
port will be described in Sec. |III A| While in principle 
exact, in practice TDDFT requires the use of approx- 
imations for the time-dependent XC functionals. The 
most popular one, the adiabatic local density approxi- 
mation (ALDA), depends only on the local and instan- 
taneous density, i.e., it does not include memory effects. 
On one hand this feature is certainly a shortcoming of the 
approximation whose consequences for time-dependent 
transport still need to be explored. On the other hand, 
in the context of multistability, this approximation allows 
one to formulate the steady-state condition in terms of 
a closed set of non-linear equations for the steady-state 
density. The solution of these equations allows for an ef- 
ficient scan of parameter space to find those parameter 
values for which multistability is possible. This steady- 
state approach with the local and adiabatic approxima- 
tion of TDDFT will be described in Sec. llHBl 

A. Real-time TDDFT for transport 

One of the technical difficulties in applying TDDFT to 
quantum transport lies in the necessity of propagating an 
infinite non-periodic system, or equivalently, a finite open 
system attached to semi-infinite leads. Here we sketch 
the technique how this can be achieved. The technical 
details of the algorithm can be found in Ref. ESI 

In a localized (site) basis, the KS Hamiltonian of the 
total system consisting of left lead, device, and right lead 
can be partitioned in a block-diagonal matrix form as 

/Hff(i) H iC \ 

H Ci H**(t) n C R , (5) 
V U RC H*g(i)/ 

where H0^(i) is the Hamiltonian of the isolated device 
and H£f (t) = H£f + Wf s (t) is the Hamiltonian of 
the isolated lead a. The time-dependent potential in the 
leads has the simple form W* s (t) = W a (t)l a where 
1 Q is the identity matrix for lead a = L,R. Finally, 
Hcq and VLac describe the coupling between lead a and 
the device. Here and in the following we use boldface 
notation to indicate matrices in one-electron labels. 

Using downfolding techniques one can derive the equa- 
tion of motion for the fc-th KS single-particle orbital pro- 
jected onto the central region, ipk,c(t), which reads 

[id t - Hgg(t)] Vfc,c (t) = / dt E* (t, t)V*,c(t) 

Jo 

+ X) H CagS,(*,0)^fc,a(0). (6) 
a 

Here, ipk,a(0) is the projection of the fc-th KS orbital on 
lead a at the initial time to = 0, g„ a (i, t') is the retarded 
Green function of the isolated lead a and the retarded 



embedding self-energy is defined as 

S&(i,f)= £ H CQg l(M')H Q c. (7) 

a=L,R 

The time-dependent density at site j of region C at zero 
temperature can be written as 

occ 

n j (t) = 2j2\^k,c(J,t)\ 2 , (8) 

k 

where the sum runs over the occupied KS orbitals and 
the prefactor is due to spin degeneracy. 

In our model the KS Hamiltonian matrix H^<£(i) has 
a tridiagonal form where the only non-vanishing entries 
are the off-diagonal matrix elements [H^?(t)]j j+ \ = 
[ H cc(*)]j+U = -Vo with j = 1, • • ■ ,N C - 1, and the 
diagonal matrix elements 

= sf(t)+v H [n 3 (t)}+v xc [n}(j,t), (9) 

with j = 1, ••• ,Nc- The second term on the r.h.s of 
Eq. Q is the Hartree potential and the third term is 
the XC potential of TDDFT for model systemsPE3 Al- 
though in our model the interaction is restricted to the 
device region only, the exact XC potential has contribu- 
tions in the adjacent lead regions as wel l and will rigor- 
ously vanish only deep inside the leads Therefore, 
already at this point we make an approximation by re- 
stricting the XC potential to the device region only. 

Of course, the exact form of i^cIt 1 ] is unknown and in 
practice one has to resort to approximations. For lattice 
systems, a local density approximation (LDA) based on 
the Bethe ansatz solution of the unif orm o ne-dimensional 
Hubbard model has been suggeste d 41 * 42 ! anc j a parame- 
terization of this Bethe ansatz LDA (BALDA) has been 
presented in Ref.|4TJ The adiabatic version 4 ^ of this func- 
tional (ABALDA) makes u xc M local in both space and 
time, i.e., v xc [n](j,t) = v xc (rij(t}). The original BALDA 
was designed for the uniform Hubbard model with Hub- 
bard interaction U and nearest-neighbor hopping V ev- 
erywhere. A modification of this functional for the case 
of a single interacting impurity site (Nc = 1) connected 
by a hopping matrix element Vu n k to the leads with hop- 
ping V has been suggested in Ref. 05l In this work we 
use this modified functional both for the case of a single 
impuritySD and also for Nc > 1 where in the latter case 
we impose the restriction that the hopping Vc between 
the sites in the central region is equal to the hopping Vu n k 
from the chain to the leads. A particularly interesting 
property of the BALDA is its discontinuity at integer val- 
ues of the occupation numberP This discontinuity has 
a fundamental impact for time-dependent transport in 
the Coulomb blockade regime and may prevent a biased 
system from evolving towards a steady stateP^ 

Finally, we note that for our Hubbard-like form of the 
interaction, where each electron interacts only with elec- 
trons of opposite spin on the same site, also the HF poten- 
tial becomes a local potential depending only on the local 
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density. In our numerical studies we will also present re- 
sults obtained within the HF approximation. 



B. Adiabatic approximation: Steady-state 
condition for the density 

Following the time evolution of the system as it is 
driven out of equilibrium by applying a bias in the leads, 
tells us if and how the system attains a steady state in 
the long-time limit. However, without doing the actual 
time propagation, we can also assume that the system 
approaches a steady state and study the consequences 
of this assumption. In the context of multistability in 
TDDFT, a particularly useful consequence is that, within 
the local and adiabatic approximation, it is possible to 
derive a self-consistency condition for the steady-state 
density. 

Using non-equilibrium Green function techniques, the 
steady-state density fij :— linit-yoo rij(t), can be obtained 
from the lesser Green function projected onto the central 
region, G cc , as 



2 / £ [ G -<">] 



j j 



(10) 



The lesser Green function can, in turn, be obtained from 

G< C M = Gg c H£< m HG£ c H, (11) 

where £< m (u;) is the lesser embedding self-energy while 
Gqc = Pec] i s ^e retarded KS Green function which, 
at the steady state, reads 



Gg c ( W )=( W -H££(n)-E£.M 



In this formula 



H 



cc {h) := Urn Hgg(t), 



(12) 



(13) 



is the asymptotic value of the KS Hamiltonian in the 
central region. Note that H^S depends on n only in the 
local and adiabatic approximation. 

For the ID tight-binding leads the retarded embedding 
self-energy has the structure 



(14) 

and therefore the self-consistency condition ( 10 ) becomes 



(/LMr L (u;)|[Gg c (u;)] 



l.J 



+/K(w)r fl (w)|[Gg c (w)] JVeJ f) ) (15) 

with the shifted Fermi function f a (uj) = f(u) — W a ) and 
W a := limj-^oo W a (t) . Since the KS Green function of 



Eq. (12) depends on n only, Eq. (15) is a set of coupled 
nonlinear equations for the steady-state density hj, j = 
1, . . . , Nc, at the sites of the device region. Due to the 



nonlinearity of these equations, more than one fixed point 
solutions may exist. Therefore, one can use Eq. (15) 



to scan the parameter space for possible multiple steady 
states. Once these values are identified we can investigate 
if and how the different steady states are attained by time 
propagation. 

We close this Section by recalling a few basic proper- 
ties of a fixed-point(FP) solution 46 which we will use in 
the following to interpret and understand our numerical 
results. Consider a system of Nc coupled equations of 
the following general form 



n 



g(n), 



(ni,n 2 , —,n No ). 



(16) 



Let the vector n be a fixed point solution of Eq. ( 16 ) and 



J the Jacobian matrix |n- A fixed point is stable if 
the modulus of all eigenvalues of J is smaller than unity, 
partially unstable if there exist at least one eigenvalue 
with modulus larger than unity (saddle point), and it is 
totally unstable if the modulus of all eigenvalues is larger 
than unity. 



IV. MANY-BODY TECHNIQUE: 
KADANOFF-BAYM EQUATIONS 

An alternative approach to TDDFT for transp ort is 
the time-dependent MBPT formulation of transport! 20 ! 21 ! 
based on the time evolution of the Green function via 
the KB equationsP2tfG2 In the MBPT, the many-body 
self-energy Smb can be approximated by a selection of 
suitable Feynman diagrams, relevant for the description 
of the main scattering processes. In earlier worlP^ we 
have found that the 2B approximation is in excellent 
agreement with accurate TD-DMRG results in the regime 
of weak to intermediate interaction strength for the An- 
derson impurity modelF^ Hence the 2B approximation 
is extremely useful for benchmarking other approxima- 
tions. The main quantity of MBPT is the Keldysh Green 
function, G(z, z'), where z and z' are time coordinates 
on the Keldysh contour To describe the electron 

dynamics of the system, the Keldysh Green function is 
propagated in time according to the KB equations. Since 
we are interested in the dynamical processes occurring 
in the central region, we can again use embedding tech- 
niques to derive the equation of motion for the Green 
function Gcc projected on the central region. This equa- 
tion reads 

[td z -H cc (z)]G cc (z,z , ) = S(z,z / ) 

dz [S om (z, z) + S MB (z, z)\ G C c(z, z') , (17) 

Jc 

where, in addition to the embedding self-energy 
X cm (z, z), we also have included a many-body self-energy 
Smb(z, z)[G] . This later quantity is a functional of 
the Green function and fulfills all basic conservation 
laws ESEil We consider only the central region to be in- 
teracting and the leads are effectively noninteracting. 



Therefore the many-body self-energy has nonvanishing 
elements only for the central region, because the dia- 
gra mmat ic expansion starts and ends with an interaction 
line™ 

Equation (17) is an exact equation for Gees provided 
an exact expression for Smb is inserted. Of course, for 
practical calculations the many-body self-energy must be 
approximated. In this paper we explicitly considered the 
following conserving approximations for the self-energy; 
the HF, 2B and GW approximations. Their diagram- 
matic representations are illustrated in Fig. [T] The im- 
plementation of Eq. (171 requires a transformation into 



equations for real times, known as K B equatio ns, which 
are then solved by time propagation ) 16 I 2 ° I 21 I 5 H From the 
knowledge of the Green function any one-particle prop- 
erty of the system can be extracted. In particular, the 
time-dependent density can be obtained from the lesser 
Green function as 



-2i[G< c (t,t + )] jj 



(18) 



In the correlated case there is no simplification such as 
Eq. ^ since there are no more well defined single-particle 
states. Also the current through lead a can be exp ressed 
in terms of the Keldysh Green functions and reads^l2U 



I a (t) = 2Re^Tr 



re 



dt[G< c (t,t)^ a (t,t) 



diG% c (t,t)X<(t,t)} 



to 



Jo 



(19) 



Besides the superscripts A (advanced), R (retarded) and 
< (lesser) we also introduced the superscripts [ and ] 
to denote the components with one time argument on 
the imagi nary axis and the other on the real axis and 
vice versa. 23 49 The trace is taken over onc-clcctron in- 



dices in the central region. Eq. ( 19 ) generalizes the Meir- 
Wingreen formula^ as it includes initial many-body and 
embedding effects through the integral along the imagi- 
nary track of the contour (last term). 

For the interpretation of the numerical results we found 
it useful to look at the TD spectral function defined ac- 
cording to 

A(T» = -ImTr J —e*^ [G> - G<] (T + ~, T - ~), 

' " (20) 
where r = t — t' is a relative time and T = (t + 1')/2 is an 
average time coordinate. In equilibrium, this function is 
independent of T and has peaks below the Fermi level at 
the electron removal energies and above the Fermi level at 
the electron addition energies. If the time-dependent ex- 
ternal field becomes constant after some switching time, 
then also the spectral function becomes independent of 
T after some transient period and has peaks at the addi- 
tion and removal energies of the non- equilibrium biased 



HF: 
2B: 

GW: 




FIG. 1: Conserving many-body approximations. 



system. 

In the HF approximation the many-body self-energy is 
frequency independent and therefore the only broaden- 
ing of the spectral peaks comes from the embedding part. 
This is also the case for the ABALDA. When going be- 
yond the mean-field level, the self-energy becomes fre- 
quency dependent and as a consequence the peaks of the 
spectral function are typically broadened. 

In the local and adiabatic approximation of TDDFT 
we have found a shortcut to determine if, for a given pa- 
rameter set, the biased system can have multiple steady 
states (see Section IIIB). It is important to note that 
a similar shortcut does not exist in MBPT. The reason 



is that in MBPT the self-consistency condition for the 
steady-state density requires the knowledge of Gcc(w) 
(to compute the many-body self-energy Smb[Gcc]) and 
thus the equations cannot be closed in terms of the den- 
sity only. As a consequence, in MBPT one does not have 
an efficient method to scan the parameter space to look 
for bistability. Although this feature is somewhat in- 
convenient for our analysis, physically it is quite reason- 
able because it implies that possible dynamical XC effects 
(arising from the frequency dependence of the many-body 
self-energy) are included. 



V. RESULTS 

In this Section, we present the results of our numerical 
simulation for a certain set of parameters for which the 
self-consistency condition ( 15 ) admits multiple solutions. 



We investigate how one can switch between different sta- 
ble solutions by applying a time-dependent gate voltage. 
We also demonstrate that for the same parameter sets 
the bistability is suppressed in the correlated many-body 
approximations, e.g., 2B and GW. The analysis will be 
carried out in two types of model devices, namely the one 
and two-site Hubbard models. 



A. Single site Hubbard model 

As a first example, we study a single-site Hubbard 
model connected to semi-infinite leads with the following 
parameters: V Vmk = 0.3, W L = 1.8, W R = -1.0, U = 
2.0, e c = -0.6, e a = Ef = (half-filled leads), and the 
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FIG. 2: Spectral functions for the different steady-state so- 
lution of the HF approximation (top-left panel) and BALDA 
(top-right panel). The graphical solution of Eq. ( |15[ ) is dis- 
played in the bottom-left panel for the HF and BALDA. For 
comparison we also report the 2B and GW steady-state spec- 
tral functions in the bottom-right panel. 



inverse temperature j3 = 90. All energies are measured 
in units of the lead-hopping parameter V. In the 
biased system the band-width of the leads extends from 
e F + W a — 2V to e F + W a + 2V. With these parameters 
the self-consistent equation ( |l5| admits five (three) 
solutions within the HF(BALDA) approximation. The 
fixed points are shown in the lower left corner of Fig. 2] 
where we display the left and right hand side of Eq. ( 15 1. 
The corresponding densities for the HF are fix = 0.17, 
fi2 = 0.54, n-3 = 1.0, fi4 = 1.46 and = 1.83 while for 
the BALDA the three fixed point densities are fi\ — 0.18, 
r~i 2 = 1.00, n 3 = 1.82. 



For a single site the fixed point theorem tells us that 
a solution is stable if 1 4 s | < 1, with g being the right- 



hand side of Eq. (151. Hence, one can see from Fig. [2] 
that the fixed points fi\ , and are stable in the case 
of the HF, while in case of the BALDA the stable solu- 
tions are h\ and n^. Although the solution with density 
of unity exists for both approximations, it is stable in 
the HF approximation and unstable in the BALDA. In 
the upper panels of Fig. [2] we plot the steady-state spec- 
tral functions corresponding to the fixed points of the 
HF and BALDA. The HF peak for density ni = 0.17 
(and the BALDA peak for density ri\ = 0.18) is located 
within the right lead energy continuum, while the HF 
peak for density — 1.83 (and the BALDA peak for 
density = 1.82) is located within the left lead energy 
continuum. By contrast, the HF spectral function of the 
unstable fixed points, n 2 and 77,4, are peaked at the edges 
of left and right lead band respectively. The HF and 
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FIG. 3: Top panel: Time-dependent density in the HF ap- 
proximation (top) and ABALDA (bottom) after the sudden 
switch on of the bias voltage and a series of gate pulses as in 
Eq. (21 1. Bottom panel: Time-dependent density within 2B 



(left) and GW approximations (right) after the sudden switch 
on of the bias voltage and a gate pulse as in Eq. (21 1 with 
V g = -3,0,3. 



BALDA spectral functions of the fixed point with den- 
sity of unity are identical (the XC potential is zero in 
this case) and the peak is located exactly in the mid- 
dle of the overlapping region between the left and right 
bands. In spite of this the stability condition of this fixed 
point is completely different in the HF and BALDA case. 
Since the multistability can be most easily observed if the 
spectral peaks of the stable solutions are well separated, 
we conclude that this phenomenon is favoured when the 
energy bands have a small overlap and the system is in 
the negative differential resistance (NDR) regime. As we 
shall see below, for the correlated MBPT approximations 
the situation is very different. In the lower left panel of 
Fig. [2] we show the 2B and GW steady-state spectral 
functions, as obtained from the propagation of the KB 
equations. The spectral weight is spread over the whole 
lead energy range and beyond. Consequently, the height 
of the spectral function is also much smaller. The con- 
siderable broadening is due to an increased quasi-particle 
scattering in the out of equilibrium system as already ob- 
served in Ref. [2"T1 

Let us now study how to switch between different sta- 
ble steady-state densities using ultrafast time-dependent 
driving fields. We start from the initially unbiased equi- 
librium system with ground-state density hq = 0.69 
(n = 0.82) for HF (BALDA). In Fig. [3] we show the 
time evolution of the density at the interacting site for 
different approximations after the sudden switch on of 
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the bias voltage W L = 1.8 and W R = -1.0. In the HF 
approximation we observe that after some transient time 
the density approaches the value 1. The behavior of the 
ABALDA density is very different, in agreement with 
the fact that the solution with density of unity is unsta- 
ble and hence cannot be reached by time evolution. At 
the steady state the ABALDA density equals the lowest 
value hi. 

To switch to the other stable solutions in real time we 
applied a time-dependent gate pulse on the Hubbard site. 
In this work we have used an exponentially decaying gate 
voltage of the form 



Vg(t) = 



V a e-* 

-V g e-^- T ^ 

V g e-^- 2T ^ 



if < t < Tg 
if T g < t < 2T g 
if t > 2T n 



(21) 



In Fig. [3] we show that in the HF case, the state with 
the lowest density hi can be obtained (in addition to 
applying a sudden bias in the leads) by switching on a 
pulse with amplitude V g = —3.0, decay rate 7 = 0.2 and 
T g = 00. The state with highest density h 5 = 1.83 can 
be obtained in a similar fashion but now applying a gate 
with positive amplitude V g = 3.0. Thus, by changing 
the amplitude of the gate voltage we can switch between 
stable steady-state solutions. For instance, with a first 
pulse of positive amplitude and T g = 50 ^> I/7 the sys- 
tem reaches the state with h§. At the time T g we apply 
a second pulse but with negative amplitude. The density 
shows a transient behavior after which it approaches the 
value hi. If we now apply a third pulse of positive am- 
plitude at time 2T g the density goes back to the initial 
value n 5 . This is nicely illustrated in Fig. [3] In Fig. [4] we 
show the non-equilibrium HF spectral function A(T, ui) of 
Eq. (20) for a double switch with V g = —3.0, 7 = 0.2 and 



T g = 50. The figure enlightens an interesting aspect re- 
garding the transition from one steady-state to another. 
The density rises from the lowest hi to the highest h§ 
lingering for a while in the middle stable solution hs . 

Going back to Fig. [3] we see that also in the ABALDA 
the state with highest density n.3 = 1.83 is reached by 
applying a gate pulse with V g — 3.0 and 7 = 0.2 (in 
addition to a sudden bias in the leads) . If the amplitude 
is negative instead (V g — —3.0) the density increases first 
but eventually drops down and goes back to its initial 
value hi. Like in the HF case we can switch back and 
forth between stable solutions by changing the sign of 
V g . Not unexpectedly, however, the decay time r 7 ~ 
I/7 cannot be arbitrarily short. If r 7 is too short the 
system does not have time to accumulate or lose enough 
density to change the self-consistent potential and after 
some transient it falls back to the previous steady-state 
value. This is clearly shown in Fig. [3] for the amplitude 
V g = ±3.0, T g = 50 and a faster decay rate 7 = 0.6. 

Intuitively one would expect that by increasing (de- 
creasing) the on-site energy of the Hubbard site the den- 
sity decreases (increases). However, the highest (lowest) 
stable steady-state density is obtained with a positive 




FIG. 4: Non-equilibrium spectral function for a gate pulse 
V g — —3, 7 = 0.2 which brings the density to h\ first, fol- 
lowed by a second identical gate pulse gate but with opposite 
amplitude which brings the density to hs. The intermediate 
transition to the n.3 stable solution is clearly visible. 



(negative) gate. This is due to the fact that in our case 
the on-site energy of the Hubbard site lies below the en- 
ergy band of the left lead. By applying a positive gate 
a finite hybridization occurs, leading to the migration of 
extra charge from the left lead to the Hubbard site. A 
similar argument explains the reduction of the density on 
the impurity site when a negative gate is turned on. 

In the lower panels of Fig. [3] we plot the densities ob- 
tained within the 2B and GW self-energy approxima- 
tions. We applied the bias voltage and a gate pulse of 
the form V g (t) = V g e~^ for t > with V g = 0, ±3. 
In all cases only one steady state emerges at the end of 
the propagation with a density of about 1.0. It is worth 
observing that the 2B and GW steady-state values of 
the densities are close to each other, indicating that the 
single-bubble diagram, common to both approximations, 
is the dominant term of the perturbative expansion in 
this casePS 

By time propagation we have shown that the three HF 
densities, hi, h% and 77.5 , and the two ABALDA densities, 
hi and 713, are stable in a slightly different sense than 
that of the fixed point theorem. The fixed point theorem 
does not contain any information on the actual dynam- 
ics. Similarly, the HF solutions 77,2 and hi as well as the 
ABALDA solution n 2 are unstable in the sense that there 
exist no external perturbation to drive the system toward 
them. Thus, the fixed point theorem provides us with a 
good criterion to establish whether a given steady-state 
can be reached or not. This criterion is certainly rigorous 
in the limit of adiabatic switchings but, as we just found, 
its validity extends well beyond the adiabatic regime. 

The time-dependent currents for the various approxi- 
mations are shown in Fig. [5] Corresponding to the three 
stable HF steady-state densities there exist only two dis- 
tinguishable values for the current I nit) at the interface 
between the Hubbard site and the right lead. The lower 
value corresponds to the solutions hi and n.5, while the 
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FIG. 5: Time-dependent current after the sudden switch on 
of the bias voltage and of a gate pulse as in Eq. ( |21[ ) in the HF 
(top-upper panel), ABALDA (top-lower panel), 2B (bottom- 
left panel) and GW (bottom-right panel) approximations. In 
the HF and ABALDA case a time-dependent switch between 
two different steady-states is shown. 



higher value corresponds to the solution n 3 . The exis- 
tence of only two solutions for the current is the conse- 
quence of an approximate particle-hole symmetry of the 
self-consistent equation (15 1, i.e., n§ ~ 1 — hi. One par- 



ticular appealing feature of the HF currents is the large 
difference between the two steady-state values, a highly 
desirable property for designing nanoscale diodes. 

The particle-hole symmetry holds also for the 
ABALDA and therefore the steady-state currents corre- 
sponding to the two stable solutions are almost indistin- 
guishable. Finally, the 2B and GW steady-state values 
of the currents, approach the same value independent of 
the gate voltage, in agreement with the existence of a 
unique steady state. 

Increasing the bias the number of HF fixed point solu- 
tions reduces to three of which only two are stable. Also 
by increasing the interaction strength U the number of 
stable solutions reduces to two because a small amount 
of density causes a considerable change in the effective 
potential. Consequently, the middle solution becomes 
unstable. 



B. Two site Hubbard model 

In this Section we consider the case of two interact- 
ing sites (Nc = 2) connected to two semi-infinite, non- 
interacting tight-binding leads. We choose the following 
parameters: V hnk = 0.4, W L = 2.2, W R = -1.2, U = 2.0, 
V G = Vi.2 = 0.4, e a = e F = 0, ef = = -0.6 and /3 = 



Upper panel: Graphical solution of the integral in 
Lower panel: Spectral functions for the HF approxi- 
mation with Hubbard interactions corresponding to the seven 
different steady-state solutions for the density. 



Eq. (151 



FP 


ni 


n 2 


FP 


rai 


n 2 


1 


0.094 


0.144 


5 


1.867 


1.862 


2 


0.150 


1.146 


6 


1.129 


0.546 


3 


0.124 


1.860 


7 


1.794 


0.226 


4 


1.098 


1.821 









TABLE I: Fixed point (FP) solutions of Eq. gHJ for the 
steady-state densities of two interacting Hubbard sites con- 
nected to two biased, non-interacting leads in the HF approx- 
imation (see upper panel of Fig. IHJ. The parameters are: Viink 
= 0.4, W L = 2.2, W R = -1.2, U = 2.0, V 1>2 = 0.4, e a = e F = 0, 
and ef = e% = -0.6. 



90. The leads are half-filled and the lead bands have an 
energy range between e a + W a — 2V and e Q + W a + 2V. 
Within the HF approximation, the steady-state con- 



dition ( 15 ) then has seven solutions which are shown in 
the upper panel of Fig. [6j The black curve is obtained 
by finding the root of the equation n-x — gi(n\, n%) = at 
fixed ni where gi (n\ ^n-i) is the right hand side of Eq. ( 15 ) 
with j = 2. The red curve is obtained in an analogous 
way by exchanging 1 <H> 2. Hence the intersections of the 
curves are the fixed points. The numerical values for the 
steady-state densities at the two Hubbard sites for the 
seven fixed points are given in Table [TJ 

In the lower panel of Fig. [6] we show the spectral func- 
tions corresponding to the seven different fixed points 
(FP's). The spectral function for FP 1 is located mostly 
in an energy range within the energy band of the right 
lead, while the one for FP 5 has most of its weight in 
the energy range of the left lead. In contrast, FP's 6 and 
7 have considerable weight in the energy bands of both 
leads. The spectral functions corresponding to FP's 2, 
3, and 4 have much narrower peaks than the spectral 
functions for the other fixed points. 

According to the fixed point theorem only the FP's 1, 
3, 5, and 7 are stable and we expect them to be accessi- 
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FIG. 7: Densities and currents for the HF approximation in 
case of short range (Hubbard) and long range (Coulomb) in- 
teractions. A switch between different steady-states by ap- 
plying exponentially decaying gate of Eq. ((21 1) are shown. 



ble by time propagation. In the upper left panel of Fig. [7] 
we show the time evolution of the total density on the 
two dots, n tot (t) = ni(t) + n 2 (t), in the HF approxima- 
tion for a sudden switch-on of the bias and several gate 
voltages starting from the equilibrium state with ground 
state density n\ — n 2 — 0.83. The steady state corre- 
sponding to FP 1 is obtained by applying only the bias 
(no gate) . In the case where we apply, in addition to the 
bias, a decaying gate voltage of the form (21 1 to both sites 
with V 9 ! = V g 2 = V g = 3.0 and decay rate 7 = 0.2, the 
total density increases and after some transient evolves 
towards the steady state corresponding to FP 5. In this 
case, lifting the on-site energy due to the gate voltage 
allows for extra charge to accumulate at the interacting 
sites such that the high-density steady-state solution can 
be achieved. In contrast, the solution of FP 7 can be ob- 
tained by applying the decaying gate voltage to the first 
(left) site only with amplitude V g \ — 3.0, V g 2 = and 
7 = 0.2. 

Surprisingly, applying a similar asymmetric gate volt- 
age but with a smaller amplitude (V 9 .i = 1.0, V g .2 = 0), 
leads to a very different long time behavior. In this case 
the system does not evolve towards a steady state after 
the transients, instead we observe an oscillatory time- 
dependent density. In the long-time limit, the time- 
dependent total density (purple curve in the upper left 
panel of Fig. [7| oscillates with an amplitude of the order 
of 10~ 3 , around 1.96. This value corresponds approxi- 
mately to the total steady-state density of FP 3 of Fig. [6] 

Despite this apparent similarity, the nature of these 
solutions is very different. While for the steady state of 
FP 7 the charge is mostly located at the first site, in the 
case of FP 3 the density on the first site is smaller than 




FIG. 8: Densities of the first and second site for the HF ap- 
proximation in case of Hubbard interactions corresponding to 
the middle curves in the upper left panel of Fig. [7J 



on the second one (see Table [T]). The different nature of 
these two cases becomes even more obvious when looking 
at the time evolution of the density at the two interact- 
ing sites separately (see Fig. [8]). While in the first case 
{Vg,i = 3.0, V 9: 2 = 0) the steady state is attained quite 
rapidly, in the second case (V g _i — 1.0, V g _2 — 0) we see 
non-decaying density oscillations at the individual sites 
with rather large amplitudes. In the long-time limit the 
density oscillates thus inducing an oscillating KS poten- 
tial. The persistence of these oscillations means that the 
density solves the Floquet system of equations in which 
the harmonics of the potential depend on the density it- 
self. At first sight one might be reminded of non-decaying 
density and current oscillations which can appear for non- 
interacting systems when the biased system possesses two 
or more bound states. However, here we work in the HF 
approximation and therefore the analysis of Refs. I51H561 
needs to be modified, see below. 

Some insight into the nature of these oscillations can 
be gained from the simple model of an isolated Hubbard 
dimer. In the HF approximation, the equation of motion 
for the electronic density matrix p of the isolated dimer 
reads 



id tP (t) = [H HF (t),p(t)} ) 
where the HF Hamiltonian is given by 



H HF (t) 



ei + Um/2 Vi.2 

e 2 + Un 2 /2 



(22) 



(23) 



Under the simplifying assumption e\ — £2 one can then 
derive the equation of motion for the quantity 5n(t) = 
ni(t) — ri2(t) = pn(t) — P22(t) which reads 

772 

dn + (W? 2 - UD)Sn + — (Sn) 3 = 0, (24) 

8 

where the constant D is related to the initial condition 



of Eq. ( 23 ) and can be defined through the off-diagonal 
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FIG. 9: Oscillation frequency and amplitude of the density 
oscillations found in HF for certain gate switchings of the 
Hubbard dimer connected to biased leads as function of the 
hopping between the Hubbard sites. For comparison, oscil- 
lation frequencies and amplitudes are given for the isolated 
Hubbard dimer in HF approximation described by Eq. (24 1. 



matrix elements of the density matrix as 

D = Vi, 2 (pi :2 (0) + ,92,1(0)) + ?(£n(0)) 



(25) 



We note that Eq. (24 1 is the equation of motion of a 
classical, anharmonic oscillator and therefore supports 
oscillating solutions. We now check if the model of the 
isolated Hubbard dimer has anything to do with the os- 
cillations seen in our tra nsport setup. To this end, we 

calculate D from Eq. (24), i.e., D = ^ + ^ + f(<5™) 2 , 
where the densities and their time-derivatives are taken 
from the transport calculation after the transients have 



died out. As the Eq. (24) is an approximation for the 
connected Hubbard dimer, D is not constant in time. 
Hence in order to compare the oscillation amplitudes and 
frequencies from the transport simulations with those re- 
sulting from Eq. p4| we averaged D over an oscillation 
period. In Fig. [9] we show the dependence of oscillation 
frequency and amplitude for different switchings of the 
gate {V g<1 (t) = V exp(-yt), V g , 2 (t) = 0, V = 1) as 
function of the hopping Vi,2 between the two sites of the 
Hubbard dimer connected to biased leads and compare to 
the corresponding solutions of Eq. ( 24 ) . We see that both 



frequency and amplitude of the isolated and connected 
dimer behave qualitatively quite similar as function of 
the intersite hopping and we conclude that the model 
of the isolated dimer certainly captures the physics be- 
hind these oscillations. We also would like to point out 
that the regions of parameter space where the oscillations 
are found appears to be quite small. For most parame- 
ters the system actually does evolve towards one of the 
steady states of Table [TJ 

The occurrence of self-induced persistent oscillations 
in the HF mean field theory is favoured by the short- 
range nature of the Hubbard interaction. In fact, we also 
have studied a modified version of our model where we 
replaced the last term of Eq. ^ by a more long-range, 
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n 2 


FP 
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ri2 
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0.147 
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1.585 


0.674 
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0.632 


1.658 
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1.624 


0.250 
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1.506 
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TABLE II: Fixed point (FP) solutions of Eq. (jT5j) for the 
steady-state densities of two interacting Hubbard sites con- 
nected to two biased, non-interacting lead s in the BALDA 
approximation (see upper panel of Fig. 111. The parameters 
are: Vn nk = 0.4, Wl = 2.2, Wn = -1.2, U = 2.0, Vi, 2 = 0.4, 
e a = e F = 0, and ef = -0.04, e% = 0.2. 



Coulomb-like interaction i J^t 



t it 



d) al dj a 'd l(y with 



U 

u 

, 2|i-j| 



* = 3 

i ^ 3 



(26) 



In this case we have not found any oscillating solutions 
in the long-time limit. We have found two stable steady- 
state solutions accessible by time propagation. The first 
steady-state solution has densities ri\ = 1.06, n 2 = 1.09, 
the second one has n\ — 0.11, n 2 = 0.13. The spectral 
functions corresponding to these solutions (see Fig. 10) 



arc localized around the Fcrmi-lcvel of the left or right 
lead respectively. The inclusion of the long range interac- 
tion destroys the states where the first peak is localized 
on the right lead energy band and the second peak is 
localized to the left lead energy band. Because the mag- 
nitude of the interaction felt by the electron on the site 
is now higher the density on the sites is decreased and 
the solution corresponding to the highest density is at 
half-filling. Also in this case we are able to switch be- 
tween the two steady-states. The currents corresponding 
to these two solutions for the density have almost the 
same magnitude. 



In Fig. 10 we show the time-dependent densities and 
currents for the 2B approximation. Again within the cor- 
related approximations we find only one solution for the 
density and current. In the lower panels of Fig. 10 we 
show the spectral functions for the 2B and GW approx- 
imations compared to the spectral functions of the HF 
approximation. The 2B and GW spectral functions are 
qualitatively quite different from those of the HF approx- 
imation. Instead of the two peak structure of HF approx- 
imation, with 2B and GW approximations we have one 
very broad peak with much lower maximum. 

We also studied the possibility of multiple steady- 
states for the same model within the BALDA. Using the 
same parameters as above, the BALDA has multiple so- 
lutions. However, at least one fixed point has a density 
on one of the dots very close to unity, exactly where the 
BALDA potential is discontinuous. For a single interact- 
ing dot, this discontinuity has been shown to be closely 
related to the Coulomb blockade phenomenon.^ For the 
purposes of the present work we avoid the regime of in- 
teger occupancy in an ABALDA treatment by changing 
the on-site energies of the interacting sites in an asym- 
metric way such that ef = —0.04 and e 2 = 0.2. With 
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FIG. 10: Densities and currents for the 2B approximation 
short range (Hubbard) and long range (Coulomb) interactions 
are used. Lower panel: Spectral functions for the different 
approximations at the end of the time propagation. 
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FIG. 11: Upper panel: Graphical solution of the integral in 
Eq. (15 I. Lower panel: Spectral functions for the BALDA 



with Hubbard interactions corresponding to the five different 
steady-state solutions for the density. 
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FIG. 12: Time-dependent density and current for the 
ABALDA with different applied gates. 



at t = 0, the system approaches the third solution if no 
external gate voltage is turned on. On the other hand, 
if a gate voltage of the form (21 1 is applied only to the 

-2.0 and 7 



second site, with amplitude V 9t 2 = —2.0 and 7 = 0.2, 
the system attains a steady-state with a density corre- 
sponding to FP 5. As before switching between these 
two steady-state density is possible by changing the sign 
of the applied gate. A similar gate voltage applied only to 
the second site and smaller amplitude (V g = —1.0) leads 
to an oscillatory time-dependent density, whose average 
total density is close to the one of FP 1 . 

The frequency of the time-dependent density is u) = 
2.24 which is close to the energy difference (lu — 2.34) 
between the peaks of spectral function. Hence, these os- 
cillations are due to the existence of bound-states in the 
biased interacting system. One possible way to explain 
the role of bound-states in the biased KS Hamiltonian is 
that in the long-time limit the KS potentials are time de- 
pendent (with the bound state eigenenergy differences as 
prominent frequencies) leading again to time-dependent 
currents and densities (by virtue of Floquet's theorem). 
This is evidently achieved in adiabatic approximations 
with the XC potential depending only on the local den- 
sity such as the ABALDA and HF approximation. 



VI. CONCLUSIONS 



these modifications, the two coupled equations given by 
Eq. (15), are solved simultaneously, yielding five fixed- 



pointsfsee Fig. 11 and Table [HJ. Among these five fixed 
points, FP 1, FP 3, and FP 5 are stable, the other two 
unstable. 

The spectral functions corresponding to FP 3 and FP 
5 have two well separated smooth peaks, while the one 
corresponding to FP 1 has two sharp peaks, the first one 
located at wi = —0.168 outside the energy range of the 
left lead, the second one at uj-i = 2.16 outside the energy 
range of the right lead. 

Again, the stable solutions are accessible by time prop- 
agation. Upon application of a sudden bias in the leads 



In this paper we have investigated by means of real- 
time propagation within MBPT and TDDFT, the exis- 
tence of multiple steady states in single and double in- 
teracting quantum dot systems connected to semi-infinite 
leads. Within the framework of MBPT one can solve the 
steady-state MBPT equation without necessarily going 
through the whole time propagation. This can be done, 
by an iterative procedure like in Ref. 1571 for example. In 
this case, starting from different initial guesses for the 
Green function, bistability would manifest itself in the 
convergence to more than one self-consistent solution. 
The advantage of the KB equations over the steady-state 
MBPT equations is that during the transient regime the 
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Green function explore a finite portion of the domain of 
all possible Green functions. Thus a single time prop- 
agation is similar to iteratively solving the steady-state 
MBPT equations for a large number of initial guesses. 

In order to find the parameter region for bist ability 
we first solved the self-consistent steady-state equations 
within the HF and BALDA approximation and deter- 
mined the regime for which multiple solutions occur. We 
show that only the stable solutions are accessible by time 
propagation. Moreover, we find that by superimposing 
an exponentially decaying gate voltage pulse to the ex- 
ternal bias, it is possible to reach the various stable solu- 
tions and also to switch between them. For the same pa- 
rameters and driving fields, we then included dynamical 
XC effects by solving the Kadanoff-Baym equations with 
MBPT self-energies at the 2B and GW level of approx- 
imation. In all studied cases where adiabatic DFT and 
HF theory predict bistability dynamical XC effects de- 
stroy the phenomenon. Here we emphasize that we have 
performed 2B and GW calculations for many more pa- 
rameter sets than those for which we have shown results 
in the present work. We have found no indication for the 
existence of multiple steady states for any of these sets. 
However, due to the vastness of the parameter space, 
we cannot rule out completely the possibility of multiple 
steady states when dynamical XC effects are included. 

We wish to point out that even though ABALDA al- 
ready contains correlations it is based on two approxima- 
tions: the local and the adiabatic one. For any non-local 



but adiabatic approximation to the TDDFT function- 
als one could still derive a self-consistency condition for 
the steady-state density in the form of coupled, nonlin- 
ear equations. Because of this nonlinearity multiple solu- 
tions, i.e., multiple steady states, can be possible. There- 
fore our results suggest that it is the adiabatic approx- 
imation which permits bistability while we expect that 
the inclusion of memory effects suppresses it. We also 
conclude that bistable regimes induced by the electron- 
electron interaction only, are unlikely to be found in Hub- 
bard or extended Hubbard model nanoj unctions, and 
that other degrees of freedom, like molecular vibrations 
or nuclear coordinates, must be taken into account. 
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